// Trabalho Final de Programacao Concorrente
// Jacobi-Richardson
// Daniel Hatanaka, 7277579
// Guilherme Kalatzis, 7239426
// Thiago Tanaka, 7153014
// Prof. Dr. Julio Estrella, Estagiario PAE Luiz Henrique Nunes

//export PATH=/Developer/NVIDIA/CUDA-5.0/bin:$PATH
//export DYLD_LIBRARY_PATH=/Developer/NVIDIA/CUDA-5.0/lib:$DYLD_LIBRARY_PATH

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<cuda.h>

int J_ORDER;
int J_ROW_TEST;
int J_ITE_MAX;
float J_ERROR;

// Definindo o número de Threads por bloco
#define NTHREADS 64

// Função KERNEL do CUDA para resolver o sistema representado por uma matriz pelo método de Jacobi, cada BLOCO é responsável por calcular
// um valor de X, exemplo bloco 0 calcula calor de X0, bloco 1 o valor de X1 e assim por diante.
// Cada thread calcula o valor de um elemento do somatório do método de Jacobi
// Recebe como entrada:
//	Tamanho da matriz
//	Matriz A* = (cada elemento da matriz A) / (pelo elemento da diagonal da sua linha)
//	Vetor B* = (cada elemento de B) / (pela diagonal da matriz A)
//	Vetor xk = vetor com os valores de X a serem calculados
//	Vetor xk1 = vetor com os valores de X da iteração anterior
//	Vetor valor = vetor para fazer um DEBUG na placa de vídeo com suporte a cuda inferior a 3.0
__global__ void jacobi(int J_ORDER, float *MatrixAStar, float *arrayBStar, float *xk, float *xk1, float *valor) {
	int j = threadIdx.x, k = threadIdx.x, i;
	int tidblock = blockIdx.x;
	float aux = 0;
	__shared__ float value[NTHREADS];

	// Como não é possivel alocar um número grande de Threads dentro de um bloco então cada Thread calcula mais de um elemento
// do somatório
	while(k < J_ORDER) {
		// Uma variável armazena temporariamente a soma de um ou mais elementos do somatório
		aux += MatrixAStar[tidblock * J_ORDER + k] * xk[k];
		
		// DEBUG
		//valor[j] = value[j];
		//printf("%f\n", aux);
		
		// K recebe o valor do próximo elemento que será calculado
		k += NTHREADS;
	}	        
	
	// Valor dos elementos do somatório são guardados em um vetor
	value[j] = aux;
	
	// DEBUG
	//valor[j] = aux;
	
	// Sincroniza as Threads para que o vetor esteja sincronizado
	__syncthreads();

	// Calculando o somatório de dois em dois elmentos
	i = NTHREADS / 2;
	// Enquanto o não somar de dois em dois
        while(i != 0) {
		// A soma das posições depende do valor da Thread, a Thread zero que realiza a ultima soma e guarda na posição ZERO o valor
// do somatório
		if(j < i) {
			// Realiza a soma de duas posições
			value[j] += value[j+i];
			// Sincroniza as Threads para garantir que a soma de uma posição já foi realizada para continuar fazendo o somatório
			__syncthreads();
		}
		// A cada iteração diminui pela metade as Threads que fazem a soma
		i = i/2;
	}
        
	// A Thread ZERO calcula o valor de X
	if(j == 0) {
		xk1[tidblock] = arrayBStar[tidblock] - value[j];
		
		// DEBUG
		//valor[tidblock] = tidblock;
		//valor[tidblock] = value[j];
	}
}




int main(){
	// Ponteiro para o arquivo
	FILE *pFile;

	// Contatores
	int i;
	int j;
	int k;
	// Variárias que auxiliam no calculo pelo método de Jacobi
	float diag;
	float resultado;
	float erro;
	float maxNumerador;
	float maxDenominador;
	// Ponteiros para alocação dinamica do conteudo vindo do arquivo
	float **MatrixA;
	float *arrayB;
	float *MatrixAStar;
	float *arrayBStar;
	float **x;
	float *diff;
	// Ponteiros para alocação dinamica na memória da placa gráfica
	float *cudaMatrixAStar;
	float *cudaarrayBStar;
	float *cudaxk;
	float *cudaxk1;
	// Variáveis para calculo dos tempos
	float tempo1;
	float tempo2;
	float tempo3;
	// Variáveis para calculo dos tempos da biblioteca cuda.h
        cudaEvent_t inicial;
        cudaEvent_t tempo;
        cudaEvent_t inicialAloc;
        cudaEvent_t tempoAloc;
        cudaEvent_t inicialInsercao;
        cudaEvent_t tempoInsercao;
        // Variável de erro (para melhor debugar o programa)
        cudaError_t error;
        
	// Cria os eventos em CUDA
        error = cudaEventCreate(&inicial);
        if (error != cudaSuccess){
        	fprintf(stderr, "Failed to create INICIAL event (error code %s)!\n", cudaGetErrorString(error));
        	exit(EXIT_FAILURE);
    	}
        
	error = cudaEventCreate(&tempo);
	if (error != cudaSuccess){
        	fprintf(stderr, "Failed to create TEMPO event (error code %s)!\n", cudaGetErrorString(error));
        	exit(EXIT_FAILURE);
    	}
		
	error = cudaEventCreate(&inicialAloc);
	if (error != cudaSuccess){
        	fprintf(stderr, "Failed to create INICIAL ALOC event (error code %s)!\n", cudaGetErrorString(error));
        	exit(EXIT_FAILURE);
    	}		
		
	error = cudaEventCreate(&tempoAloc);
	if (error != cudaSuccess){
        	fprintf(stderr, "Failed to create TEMPO ALOC event (error code %s)!\n", cudaGetErrorString(error));
        	exit(EXIT_FAILURE);
    	}
    			
	error = cudaEventCreate(&inicialInsercao);
	if (error != cudaSuccess){
        	fprintf(stderr, "Failed to create INICIAL INSERCAO event (error code %s)!\n", cudaGetErrorString(error));
        	exit(EXIT_FAILURE);
    	}
    			
	error = cudaEventCreate(&tempoInsercao);
	if (error != cudaSuccess) {
        	fprintf(stderr, "Failed to create TEMPO INSERCAO event (error code %s)!\n", cudaGetErrorString(error));
        	exit(EXIT_FAILURE);
    	}		

	// Variável para armazenar o nome do arquivo
        char fileName[100];
	puts("Digite o nome do arquivo:");
	// Lê o nome do Arquivo com os dados
	scanf("%s", fileName);
	// Abre o arquivo
	pFile = fopen(fileName, "r");
	// Se não conseguiu abrir o arquivo
	if (pFile == NULL) {
		// Mensagem de erro
		fputs("Read file error\n", stderr);
		// Sai do programa
		exit(1);
	}
	else {
		// Lê os dados
		fscanf(pFile, "%d %d %f %d", &J_ORDER, &J_ROW_TEST, &J_ERROR, &J_ITE_MAX);
		// Evento para calcular o tempo para alocar as variáveis
        	cudaEventRecord( inicialAloc, 0 );
		// Aloca as variáveis na memória principal do computador
        	arrayB = (float*)calloc(J_ORDER,sizeof(float));
       	 	arrayBStar = (float*)calloc(J_ORDER,sizeof(float));
        	diff = (float*)calloc(J_ORDER,sizeof(float));
        	MatrixA = (float**)calloc(J_ORDER ,sizeof(float*));
        	MatrixAStar = (float*)calloc(J_ORDER * J_ORDER,sizeof(float));
        	for(i = 0; i < J_ORDER; i++)
                	MatrixA[i] = (float*)calloc(J_ORDER,sizeof(float));
                	
        	x = (float**)calloc((J_ITE_MAX + 2),sizeof(float*));
        	for (i = 0; i < J_ITE_MAX + 2; i++)
               		x[i] = (float*)calloc(J_ORDER,sizeof(float));

		// Aloca a variável na memória da placa gráfica e verifica se ocorreu algum erro
  		error = cudaMalloc((float**)&cudaMatrixAStar, J_ORDER * J_ORDER * sizeof(float));
    		if (error != cudaSuccess) {
        		printf("cudaMalloc MatrixAStar returned error code %d, line(%d)\n", error, __LINE__);
        		exit(EXIT_FAILURE);
    		}
		// Aloca a variável na memória da placa gráfica e verifica se ocorreu algum erro
        	error = cudaMalloc((float**)&cudaxk1, J_ORDER * sizeof(float));
        	if (error != cudaSuccess) {
        		printf("cudaMalloc cudaxk1 returned error code %d, line(%d)\n", error, __LINE__);
        		exit(EXIT_FAILURE);
    		}
		// Aloca a variável na memória da placa gráfica e verifica se ocorreu algum erro
        	error = cudaMalloc((float**)&cudaarrayBStar, J_ORDER * sizeof(float));
        	if (error != cudaSuccess) {
        		printf("cudaMalloc cudaarrayBStar returned error code %d, line(%d)\n", error, __LINE__);
        		exit(EXIT_FAILURE);
    		}
		// Aloca a variável na memória da placa gráfica e verifica se ocorreu algum erro
        	error = cudaMalloc((float**)&cudaxk, J_ORDER * sizeof(float));
        	if (error != cudaSuccess){
        		printf("cudaMalloc cudaxk returned error code %d, line(%d)\n", error, __LINE__);
        		exit(EXIT_FAILURE);
    		}
        	
		// Evento para calcular tempo de alocação
        	cudaEventRecord(tempoAloc, 0);
		cudaEventSynchronize(tempoAloc);
		
		// Calculo do tempo de alocação
		cudaEventElapsedTime(&tempo1, inicialAloc, tempoAloc);

		// Imprime o tempo de alocação
		printf("----------------------------------------------------------------\n");
        	printf("Tempo de alocacao: %.10f Miliseg\n", tempo1);
        	printf("----------------------------------------------------------------\n");
			
		// Evento para calcular o tempo de inserção dos dados
		cudaEventRecord(inicialInsercao, 0);

		// Insere os dados da Matriz A
		for(i = 0; i < J_ORDER; i++)
			for(j = 0; j < J_ORDER; j++)
				fscanf(pFile,"%f", &MatrixA[i][j]);
		
		// Insere os dados do Vetor B
		for(i = 0; i < J_ORDER; i++)
			fscanf(pFile,"%f", &arrayB[i]);
			
		// Evento para calcular o tempo de inserção dos dados
		cudaEventRecord(tempoInsercao, 0);
		cudaEventSynchronize(tempoInsercao);

		// Calculo do tempo de inserção
		cudaEventElapsedTime(&tempo2, inicialInsercao, tempoInsercao);
		
		// Imprime o tempo de inserção
		printf("----------------------------------------------------------------\n");
        	printf("Tempo de insercao de dados: %.10f Miliseg\n", tempo2);
        	printf("----------------------------------------------------------------\n");
		
		// Evento para calcular o tempo de execução
		cudaEventRecord(inicial, 0);

		// Calcula o Vetor B*
        	for(i = 0; i < J_ORDER; i++) {
                	diag = MatrixA[i][i];
                	arrayBStar[i] = arrayB[i] / diag;
       		 }
	
		// Calcula a Matriz A*
        	for(i = 0; i < J_ORDER; i++){
                	diag = MatrixA[i][i];
                	k = 0;
                	for(j = 0; j < J_ORDER; j++)
                        	if(i != j)
					MatrixAStar[i * J_ORDER + j] = MatrixA[i][j]/diag;
				else
					MatrixAStar[i * J_ORDER + j] = 0;
        	}

		// Copia a Matriz A* da memória principal para a memória da placa gráfica
        	cudaMemcpy(cudaMatrixAStar, MatrixAStar, J_ORDER * J_ORDER * sizeof(float),cudaMemcpyHostToDevice);
		// Copia o Vetor B* da memória principal para a memória da placa gráfica
        	cudaMemcpy(cudaarrayBStar, arrayBStar, J_ORDER * sizeof(float),cudaMemcpyHostToDevice);
		
		// Calcula o X inicial
        	for(i = 0; i < J_ORDER; i++)
                	x[0][i] = 0;
         
		// Define o número de Blocos por Grid de acordo com o tamanho da ordem da Matriz
        	dim3 grid(J_ORDER, 1, 1);
		// Define o número de Threads por Bloco de acordo com um número definido no começo do programa
        	dim3 thread(NTHREADS, 1, 1);
        	
		// DEBUG
       	 	float *valor1, *valor2;
        	valor1 = (float *)malloc(sizeof(float) * J_ORDER);
        	cudaMalloc((float**)&valor2, sizeof(float) * J_ORDER);

		// Calcula os valores de X até um número máximo de iterações
        	for(k = 0; k < J_ITE_MAX; k++) {
			// Copia os valores do X da iteração anterior para memória da placa gráfica
			cudaMemcpy(cudaxk, x[k], J_ORDER * sizeof(float),cudaMemcpyHostToDevice);
			// Calcula o valor de X
	                jacobi<<<grid,thread>>>(J_ORDER,cudaMatrixAStar,cudaarrayBStar,cudaxk,cudaxk1, valor2);
			// Copia os valores de X da iteração atual da memória da placa gráfica para a memória principal
	                cudaMemcpy(x[k+1], cudaxk1, J_ORDER * sizeof(float),cudaMemcpyDeviceToHost);
			
			// DEBUG
			cudaMemcpy(valor1, valor2, sizeof(float) * J_ORDER,cudaMemcpyDeviceToHost);
			
			// DEBUG
			/*
			if(k == 1) {
				for(i = 0; i < NTHREADS; i++){
					printf("%d %f\n", i, valor1[i]);
				}
				printf("\n");
			}
			*/

			// Calcula a diferença dos valores do X da iteração atual com a anterior
        	        for(i = 0; i < J_ORDER; i++) {
                       		// Garantir que a diferença tenha valor positivo
				if((x[k+1][i] - x[k][i]) > 0){
					diff[i] = (x[k+1][i] - x[k][i]);
                        	}
				else{
                        		diff[i] = (x[k+1][i] - x[k][i]) * (-1);
                        	}
               	 	}

			// Numerador para calculo do Erro com valor positivo
			maxNumerador = fabs(diff[0]);
			// Denominador para calculo do Erro com valor positivo
               		maxDenominador = fabs(x[k + 1][0]);
			// Para achar a maior diferença
			for(i = 0; i < J_ORDER; i++) {
                		if(fabs(diff[i]) > maxNumerador)
                        		maxNumerador = fabs(diff[i]);
                		if(fabs(x[k+1][i]) > maxDenominador)
                        		maxDenominador = fabs(x[k+1][i]);
                	}

			// Calcula o Erro
                	erro = maxNumerador / maxDenominador;
			
			// DEBUG
			//printf("%d  erro: %f\n", k, erro);
			
			// Se o erro for menor que o erro desejado para de calcular X
 	               	if(erro <= J_ERROR)
                     		// Para a iteração
				break;
		}

		// Calcula o resultado aproximado de uma dada equação
        	for(i = 0; i < J_ORDER; i++)
                	resultado = resultado + MatrixA[J_ROW_TEST][i] * x[k+1][i];
		
		// Imprime os resultados obtidos
		printf("----------------------------------------------------------------\n");
        	printf("Iterations: %d\n", k);
        	printf("RowTest: %d => [%f] =? %f\n", J_ROW_TEST, resultado, arrayB[J_ROW_TEST]);
        	printf("----------------------------------------------------------------\n");

		// Evento para calcular o tempo de execução
		cudaEventRecord( tempo, 0 );
		cudaEventSynchronize( tempo );
		// Calcula o tempo de execução do programa
		cudaEventElapsedTime( &tempo3, inicial, tempo );

		// Imprime os resultados
		printf("----------------------------------------------------------------\n");
        	printf("Tempo de execução do Jacobi-Richardson: %.10f Miliseg\n", tempo3);
	        printf("----------------------------------------------------------------\n");
        
		// Libera a memória alocada da placa gráfica
		cudaFree(cudaMatrixAStar);
       		cudaFree(cudaarrayBStar);
       		cudaFree(cudaxk);
        	cudaFree(cudaxk1);
        
		// Destroi os eventos para calculo dos tempos
		cudaEventDestroy(inicial);
		cudaEventDestroy(tempo);
		cudaEventDestroy(inicialAloc);
		cudaEventDestroy(tempoAloc);
		cudaEventDestroy(inicialInsercao);
		cudaEventDestroy(tempoInsercao);
       
		// Libera a memória principal
		free(arrayB);
        	free(arrayBStar);
        	free(diff);
        	for (i = 0; i < J_ORDER; i++)
                	free(MatrixA[i]);
        	free(MatrixAStar);
	
		for (i = 0; i < J_ITE_MAX+2; i++)
                	free(x[i]);
        	free(x);

		// Finaliza a execução da placa gráfica
		cudaDeviceReset();
	}
	// Finaliza o progrma
	return 0;
}
